# Return calibrated EMEP4UK ammonia concentrations (2019 - 2021)
# aggregate to 3 year mean
r_nh3_conc_19_21 <- conc_nh3_r(yrs = 2019:2021) %>% app(., mean, na.rm = T)
crs(r_nh3_conc_19_21) <- "EPSG:27700"

## project to 3857
r_conc_19_21_3857 <- project(r_nh3_conc_19_21, "EPSG:3857", res = 1000) 

cols <- c("#9ce4ff","#fbff87","#ffa754","#ff5754","#b36bff")
brks <- c(0,1,2,3,4,99)
labs <- c("0 - 1","> 1 - 2","> 2 - 3","> 3 - 4","> 4")
pal <- colorBin(palette = cols,
                   bins = brks,
                   domain = values(r_conc_19_21_3857),
                na.color = "transparent")

leaflet() %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addRasterImage(x = r_conc_19_21_3857, colors = pal, project = F) %>% 
addLegend(colors = cols, 
          labels = labs, values = values(r_conc_19_21_3857))

# Return site using a site ID
sf <- getSITE(ste_type = "SSSI", id = "1")

# Search for the correct site
sf <- getSITE(ste_type = "SSSI",grep_txt  = "bog")

# Return sites by country
sac_4326 <- getSITE(ste_typ = "SAC",country  = "Northern Ireland") %>% 
  st_transform(.,4326)
spa_4326 <- getSITE(ste_typ = "SPA",country  = "Northern Ireland") %>% 
  st_transform(.,4326)

leaflet() %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addPolygons(data = sac_4326, color = "red", fillOpacity = 0.3,
weight = 1) %>%
addPolygons(data = spa_4326, color = "blue", fillOpacity = 0.3,
weight = 1) %>% 
addLegend(colors = c("red","blue"), labels = c("SAC","SPA"))

# Return site using a site ID
sf <- getSAC(id = "UK0016599")

# Search for the correct site
sf <- getSAC(grep_txt  = "Bal")

# Return English SACs
sf <- getSAC(country = "England")

# Return site using a site name and add to Leaflet map
sf_4326 <- getSAC(grep_txt  = "Ballynahone Bog") %>% 
  st_transform(.,4326)

leaflet() %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addPolygons(data = sf_4326, color = "red", fillOpacity = 0.3,
weight = 1)

# Return site using a site ID
sf <- getSPA(id = "UK9020330")

# Search for the correct site
sf <- getSPA(grep_txt  = "Solen")

# Return English SPAs
sf <- getSPA(country = "England")

# Return site using a site name and add to Leaflet map
sf_4326 <- getSPA(grep_txt  = "Solent and Dorset Coast") %>% 
  st_transform(.,4326)

leaflet() %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addPolygons(data = sf_4326, color = "red", fillOpacity = 0.3,
weight = 1)

# Return site using a site ID
sf <- getSSSI(id = "1")

# Search for the correct site
sf <- getSSSI(grep_txt  = "A'")

# Return English SSSIs
sf <- getSSSI(country = "England")

# Return site using a site name and add to Leaflet map
sf_4326 <- getSSSI(grep_txt  = "A' Mhoine") %>% 
  st_transform(.,4326)

leaflet() %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addPolygons(data = sf_4326, color = "red", fillOpacity = 0.3,
weight = 1)

# Return site using a site ID
sf_bally <- getSITE(ste_type = "SAC",id = "UK0016599")

# estimate cell coverage at 100 x 100 m
r <- cll_cvr_sf(sf = sf_bally, res = 100, prcsn = 25)

# check area matches sf object
# multiply total by res^2 to go from prop of grd sq to m2
r_tot <- (global(r,"sum")$sum * 100^2)
sf_tot <- st_area(sf_bally) %>% as.numeric()

# difference in m2
r_tot - sf_tot

# improve prscn to reduce difference in area
r <- cll_cvr_sf(sf = sf_bally, res = 100, prcsn = 2)

r_tot <- (global(r,"sum")$sum * 100^2)
sf_tot <- st_area(sf_bally) %>% as.numeric()

# difference in m2
r_tot - sf_tot

plot(r)

dt <- cll_cvr_tbl(ste_type = "SAC",id = c("UK0016599","UK0030093"),prcsn = 25,res = 1000,as_dt = F)

## improve accuracy by lowering prcsn to 2m

dt <- cll_cvr_tbl(ste_type = "SAC",id = c("UK0016599","UK0030093"),prcsn = 2,res = 1000,as_dt = F)

head(dt,10)

v <- sac_description(id = "UK0012720")
paste0(v, collapse = "\n\n - ") %>% paste0("- ",.) %>% cat(., sep = "")

dt <- sac_features(id = "UK0012720")

dt

## Return latest 3 year mean (2019 - 2021)
dt <- ste_nh3_concs(ste_type = "SAC",id = "UK0016599",as_dt = T, prcsn = 1)

dt

## Return annual estimates for 2002 - 2021
dt <- ste_nh3_concs(ste_type = "SAC",id = c("UK0016599"),
                    yrs = 2002:2021, agr_yrs = F,as_dt = T, prcsn = 5)

p <- ggplot()

p <- p + geom_point(data = dt, aes(x = Year, y = mean_nh3_conc, color = SITENAME),
                    size = 2)

p <- p + geom_line(data = dt, aes(x = Year, y = mean_nh3_conc, color = SITENAME),
                    linewidth = 0.5)

p <- p + geom_errorbar(data = dt, aes(x = Year, ymin = min_nh3_conc, ymax = max_nh3_conc,
                                      color = SITENAME),
                     alpha = 0.6)


p + ylab(bquote("Ammonia Concentration"~"("*mu*g~NH[3]~m^-3*")"))

dt <- ste_nh3_concs(ste_type = "SAC",id = c("UK0016599","UK0030093","UK0012720"),
                    yrs = 2002:2021, agr_yrs = F,as_dt = T, prcsn = 5)

p <- ggplot()

p <- p + geom_point(data = dt, aes(x = Year, y = mean_nh3_conc, color = SITENAME),
                    size = 2)

p <- p + geom_line(data = dt, aes(x = Year, y = mean_nh3_conc, color = SITENAME),
                    linewidth = 0.5)



p + ylab(bquote("Ammonia Concentration"~"("*mu*g~NH[3]~m^-3*")")) 

## Return latest 3 year mean (2019 - 2021)
dt <- ste_cle_status(ste_type = "SAC",id = "UK0016599",as_dt = T, prcsn = 1)

dt

## Return annual estimates for 2002 - 2021
dt <- ste_cle_status(ste_type = "SAC",id = "UK0016599",yrs = 2002:2021, agr_yrs = F,as_dt = T, prcsn = 1)

dt

## Return source attribution for Ballynahone Bog SAC
dt <- ste_srce_attr(ste_type = "SAC",id = "UK0016599",as_dt = T, prcsn = 25)

# select top 3 sources of N dep.
dt[Dep_typ == "Low-growing semi-natural veg",] %>% 
  .[,.(kg_ha = sum(kg_ha)), by = FOOTPRINT_NAME] %>% 
  .[order(-kg_ha),] %>% 
  .[1:3,]

## Return source attribution for Ballynahone Bog SAC
p <- ggpl_srce_attr(ste_type = "SAC",
 id = "UK0016599",
 dep_typ = "Woodland",
 sep_wet_dry = F,
 sep_nox_nhx = F,
 txt_size = 10,
 prcsn = 1)

p

## Return source attribution for Ballynahone Bog SAC for wet/dry, nox/nhx
p <- ggpl_srce_attr(ste_type = "SAC",
 id = "UK0016599",
 dep_typ = "Woodland",
 sep_wet_dry = T,
 sep_nox_nhx = T,
 txt_size = 10,
 prcsn = 1)

p